def load_netcdf_grid(path = "Argo_South_60", resolution = 1, target_var = "temp",
with_depth = False, filter_month = None, filter_season = None,
filter_start = datetime(2000,1,1), filter_end = datetime(2021,1,1), mean = None):
if with_depth:
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
gridcounts = np.zeros(shape=[15 * resolution, 360 * resolution, 200])
else:
grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
files = glob.glob("data/{}/**/*.nc".format(path), recursive=True)
#1 glob.glob - returns a list of all the files in a given folder with the extension 'nc'
#2 recursive=true - keeps on looking for all the files till the end, or the last file that exist
# creating lists
filename2 = []
nbprof2 = []
filename2_maskout = []
## End of "creating lists"
for f in tqdm(files[:]):
# tqdm - make your loops show a progress meter
d = Dataset(f)
#1 DATASET - Creating/Opening/Closing a netCDF file.
# This is the method used to open an existing netCDF file.
lat = d.variables["LATITUDE"][:]
#mask = (lat < -60) & (lat > -74.75)
mask = (lat < -60) & (lat > -74.5)
# Counting number of profiles per NetCDF files
nbprof = [len(lat)] # number of profiles per file
nbprof2.extend(nbprof)
# filename = filename
filename = ["Filename {}".format(f)]
filename2.extend(filename)
#Create a DataFrame: filename + number of profiles per file
profiles_per_files = {'Filename': filename2, 'Number of profiles per file': nbprof2}
new_df = pd.DataFrame(data=profiles_per_files)
#Save "new_df" as csv file:
new_df.to_csv('filename_nbprof.csv', index=False)
## End of "Counting number of profiles per NetCDF files"
juld = d.variables["JULD"][:]
units = d.variables["JULD"].getncattr('units')
dates = num2date(juld, units, "standard")
datemask = (dates > filter_start) & (dates < filter_end)
mask &= datemask
if filter_month:
datemask = np.array([d.month == filter_month for d in dates])
mask &= datemask
if filter_season:
if filter_season == "Summer":
datemask = np.array([d.month >= 10 or d.month <= 3 for d in dates])
elif filter_season == "Winter":
datemask = np.array([d.month > 3 and d.month < 10 for d in dates])
mask &= datemask
if any(mask):
lat = np.round((np.abs(lat[mask]) - 60) * resolution).astype(int)
lon = d.variables["LONGITUDE"][mask]
lon = np.round((lon + 180) * resolution).astype(int)
lon[lon == 360 * resolution] = 0
if with_depth:
#pres = np.round(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
pres = np.floor(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
temp = d.variables["TEMP_ADJUSTED"][mask]
if target_var == "sal":
try:
sal = d.variables["PSAL_ADJUSTED"][mask]
except:
print("No salinity for {}".format(f))
continue
else:
pres = d.variables["PRES_ADJUSTED"][mask]
temp = d.variables["TEMP_ADJUSTED"][mask, 0]
for x in np.unique(lon):
for y in np.unique(lat):
if with_depth:
ptmask = (lon == x) & (lat == y)
for z in np.unique(pres[ptmask]):
if z>= 200 or np.ma.is_masked(z):
continue
depthmask = pres[ptmask] == z
if target_var == "temp":
values_at_pt = temp[ptmask][depthmask].compressed()
elif target_var == "sal":
values_at_pt = sal[ptmask][depthmask].compressed()
if mean is not None:
mean_at_pt = mean[y, x, z]
grid[y, x, z] = np.nansum((grid[y, x, z],np.nansum(np.square(values_at_pt - mean_at_pt))))
else:
grid[y, x, z] = np.nansum((grid[y, x, z], np.nansum(values_at_pt)))
gridcounts[y, x, z] += len(values_at_pt)
else:
ptmask = (lon == x) & (lat == y)
if target_var == "n_point":
v = np.sum(ptmask)
if not np.ma.is_masked(v):
if v != 0:
#1 added a if statement checking if number of points is greater than 0
grid[y, x] = np.nansum((grid[y, x], v))
elif target_var == "temp":
temps_at_pt = temp[ptmask]
v = np.ma.mean(temps_at_pt)
# 1 MEAN - numpy.ma.mean(self, axis=None, dtype=None, out=None, keepdims=<no value>)
# Returns the average of the array elements along given axis.
# Masked entries are ignored, and result elements which are not finite
# will be masked.
sd = np.ma.std(temps_at_pt)
# numpy.std -->> Standard deviation
if not np.ma.is_masked(v):
grid[y, x] = np.nanmean((grid[y, x], v))
elif target_var == "depth":
pres_at_pt = pres[ptmask]
if len(pres_at_pt):
v = np.ma.max(pres[ptmask])
if not np.ma.is_masked(v):
grid[y, x] = np.nanmax((grid[y, x], v))
if with_depth:
grid /= gridcounts
print(np.max(gridcounts))
print(np.nanmin(grid), np.nanmean(grid), np.nanmax(grid))
if mean is not None:
grid = np.sqrt(grid)
return grid
def load_ctd_grid(resolution = 1, target_var = "temp",
with_depth = False, filter_month = None, filter_season = None,
filter_start = datetime(2000,1,1), filter_end = datetime(2021,1,1), mean = None):
if with_depth:
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
gridcounts = np.zeros(shape=[15 * resolution, 360 * resolution, 200])
else:
grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
files = glob.glob("data/ctd/**/*.nc", recursive=True)
#1 glob.glob - returns a list of all the files in a given folder with the extension 'nc'
#2 recursive=true - keeps on looking for all the files till the end, or the last file that exist
for f in tqdm(files[:]):
# tqdm - make your loops show a progress meter
d = Dataset(f)
#1 DATASET - Creating/Opening/Closing a netCDF file.
# This is the method used to open an existing netCDF file.
lat = d.variables["lat"][:]
if lat.ndim > 0:
continue
if lat > -60 or lat < -74.5:
continue
time = d.variables["time"]
units = time.getncattr('units')
date = num2date(time[:], units, "standard")
if date < filter_start or date > filter_end:
continue
if filter_month:
if date.month != filter_month:
continue
if filter_season:
if filter_season == "Summer":
if not (date.month >= 10 or date.month <= 3):
continue
elif filter_season == "Winter":
if not (date.month > 3 and date.month < 10):
continue
lat = np.round((np.abs(lat) - 60) * resolution).astype(int)
lon = d.variables["lon"][:]
lon = np.round((lon + 180) * resolution).astype(int)
if lon == (360 * resolution):
lon = 0
x = lon
y = lat
temp = d.variables["Temperature"][:]
temp[temp == -1e10] = np.nan # Way too cold!
if with_depth:
#pres = np.round(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
pres = np.floor(d.variables["Pressure"] / 10).astype(int)
if target_var == "sal":
try:
sal = d.variables["PSAL_ADJUSTED"]
except:
print("No salinity for {}".format(f))
continue
else:
pres = d.variables["Pressure"]
temp = temp[0]
if with_depth:
for z in np.unique(pres):
if z>= 200 or np.ma.is_masked(z):
continue
depthmask = pres == z
if target_var == "temp":
values_at_pt = temp[depthmask].compressed()
elif target_var == "sal":
values_at_pt = sal[depthmask].compressed()
if mean is not None:
mean_at_pt = mean[y, x, z]
grid[y, x, z] = np.nansum((grid[y, x, z],np.nansum(np.square(values_at_pt - mean_at_pt))))
else:
grid[y, x, z] = np.nansum((grid[y, x, z], np.nansum(values_at_pt)))
gridcounts[y, x, z] += len(values_at_pt)
else:
if target_var == "n_point":
grid[y, x] = np.nansum((grid[y, x], 1))
elif target_var == "temp":
temps_at_pt = temp
v = np.ma.mean(temps_at_pt)
# 1 MEAN - numpy.ma.mean(self, axis=None, dtype=None, out=None, keepdims=<no value>)
# Returns the average of the array elements along given axis.
# Masked entries are ignored, and result elements which are not finite
# will be masked.
sd = np.ma.std(temps_at_pt)
# numpy.std -->> Standard deviation
if not np.ma.is_masked(v):
grid[y, x] = np.nanmean((grid[y, x], v))
elif target_var == "depth":
pres_at_pt = pres
if len(pres_at_pt):
v = np.ma.max(pres[ptmask])
if not np.ma.is_masked(v):
grid[y, x] = np.nanmax((grid[y, x], v))
if with_depth:
grid /= gridcounts
print(np.max(gridcounts))
print(np.nanmin(grid), np.nanmean(grid), np.nanmax(grid))
if mean is not None:
grid = np.sqrt(grid)
return grid
# This version seems a lot more inefficient for some reason
def load_netcdf_grid_new(platform = "argo", resolution = 1, target_var = "temp", filter_month = None,
filter_season = None, filter_start = datetime(2000,1,1), filter_end = datetime(2020,1,1)):
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
if platform == "argo":
data = argo
elif platform == "seal":
data = seal
elif platform == "both":
data = {}
for k in argo:
data[k] = argo[k] + seal[k]
mask = (data["lat"] < -60) & (data["lat"] > -74.5)
datemask = (data["datetime"] > filter_start) & (data["datetime"] < filter_end)
mask &= datemask
if filter_month:
datemask = np.array([d.month == filter_month for d in data["datetime"]])
mask &= datemask
if filter_season:
if filter_season == "Summer":
datemask = np.array([d.month >= 10 or d.month <= 3 for d in data["datetime"]])
elif filter_season == "Winter":
datemask = np.array([d.month > 3 and d.month < 10 for d in data["datetime"]])
mask &= datemask
if any(mask):
lat = np.round((np.abs(data["lat"][mask]) - 60) * resolution).astype(int)
lon = data["lon"][mask]
lon = np.round((lon + 180) * resolution).astype(int)
lon[lon == 360 * resolution] = 0
depth = np.floor(data["pressure"][mask] / 10).compressed().astype(int)
grid[lat, lon, depth] = np.ma.mean(data["temperature"][mask])
return grid
def plot_grid(grid, resolution = 1, title="Ocean data", cmap="nipy_spectral", cbtitle="Temperature °C", vmin=None, vmax=None):
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, -80, -5))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
#plt.title("{}".format(title, 1/resolution),fontsize = 16)
#plt.title("\n{} at {}° resolution".format(title, 1/resolution),fontsize = 16)
x = np.arange(-180, 180.1, 1/resolution)
y = np.arange(-60, -75.1, 1/-resolution)
x, y = np.meshgrid(x, y)
x, y = m(x, y)
if cbtitle=="Number of profiles":
pcm = m.pcolormesh(x, y, grid, cmap=cmap, norm=colors.LogNorm(vmin, vmax))
# #norm=colors.PowerNorm(gamma=0.25))
# print("vmin",vmin,"vmax",vmax)
#ax.xaxis.grid(True, which='major')
cb = plt.colorbar(pcm, extend='max',shrink=0.785, pad=0.05 , orientation="horizontal")
#cb.locator=ticker.LogLocator(base=10)
cb.set_label(cbtitle, labelpad= 10, rotation=0, fontsize = 14)
else:
m.pcolormesh(x, y, grid, cmap=cmap, vmin=vmin, vmax=vmax) #"colorcet" "nipy_spectral" "ocean" "gist_ncar" "cubehelix" "jet" "PuOr" "RdBu" "rainbow"
#m.pcolormesh(x, y, grid, cmap=cmap, vmin=vmin, vmax=vmax) #"colorcet" "nipy_spectral" "ocean" "gist_ncar" "cubehelix" "jet" "PuOr" "RdBu" "rainbow"
cb = plt.colorbar(shrink=0.785, pad=0.05 , orientation="horizontal")
cb.set_label(cbtitle, labelpad= 10, rotation=0, fontsize = 14)
#m.pcolormesh(x, y, grid, cmap="jet", vmin=vmin, vmax=vmax)
plt.show()
def interp_nans(grid, method='linear'):
x = np.arange(0, grid.shape[1])
y = np.arange(0, grid.shape[0])
mask = np.isnan(grid)
xx, yy = np.meshgrid(x, y)
x1 = xx[~mask]
y1 = yy[~mask]
newarr = grid[~mask]
interp = interpolate.griddata((x1, y1), newarr.ravel(), (xx, yy), method=method)
if method == "linear":
interp = interp_nans(interp, "nearest")
return interp
ctd_grid = load_ctd_grid()
